Exploring the advection-diffusion equation through the subdivision collocation method: a numerical study

The current research presents a novel technique for numerically solving the one-dimensional advection-diffusion equation. This approach utilizes subdivision scheme based collocation method to interpolate the space dimension along with the finite difference method for the time derivative. The proposed technique is examined on a variety of problems and the obtained results are presented both quantitatively in tables and visually in figures. Additionally, a comparative analysis is conducted between the numerical outcomes of the proposed technique with previously published methods to validate the correctness and accuracy of the current approach. The primary objective of this research is to investigate the application of subdivision schemes in the fields of physical sciences and engineering. Our approach involves transforming the problem into a set of algebraic equations.

The Advection Diffusion Equation (ADE) is a significant mathematical model with crucial applications in various industrial, scientific, and engineering fields.It elucidates how substances or quantities disperse within a fluid medium.Prominent applications encompass environmental modeling, fluid dynamics, heat transfer, atmospheric science, pharmacokinetics, and geological processes.In several physical phenomena, the Advection Diffusion Equation delineates the combined effects of advection (bulk flow transport) and diffusion (random spreading).In engineering, it finds utility in simulating the transfer of substances like heat, contaminants, or solutes within fluids.Advection characterizes bulk movement, while diffusion accounts for the spreading driven by concentration gradients.This equation plays a pivotal role in comprehending and predicting substance dispersion in engineering applications, such as fluid dynamics, heat transfer, and pollution dispersion.Consider one-dimensional advection diffusion transport equation of the form with the initial condition in (2) and boundary conditions as stated in (3).Here, β represents the velocity of the flow, and α is the diffusion coefficient.The concentration at the point x and at time t, denoted as u(x, t), is an unknown function.The length of the channel is denoted by L, and ψ 1 (x) , φ 1 (t) , and φ 2 (t) are given functions.The governing equation (1) becomes a heat equation if the advection speed β is zero with a positive diffusion coefficient.Also ψ 1 (x) , φ 1 (t) and φ 2 (t) are the known smooth functions.
In the past, various authors have studied mathematical models from different perspectives.For instance, Mohebbi, A. and Dehghan, M., proposed a higher-order compact solution for the advection-diffusion equation

Subdivision scheme
Subdivision schemes are efficient tools for generating smooth curves and surfaces efficiently from a discrete set of control points.This is consistent and efficient iterative algorithm to be used for modeling of curve and surfaces.Subdivision schemes can be classified into two important branches such as approximating and interpolating.In the approximating scheme, the limit curve approximates the initial polygon, and after subdivision, only the newly generated control points are part of the limit curve.In the interpolating scheme, both the control points of the original control polygon and the newly generated control points lie on the limit curve after subdivision.If Q k i are the control points of the polygon at kth level and Q k+1 i are the control points at the (K + 1) th level, then the 6-points interpolating binary subdivision scheme proposed by 25 is defined as This scheme produces C 2 -continuous curve for the range 0 < ω ≤ 3 256 , particularly we choose ω = 3 256 17 , support width (−5, 5) .Also its degree generation, reproduction and approximation order is five and six respectively.Fundamental solution of scheme (4) with basis function is given in (5) and ( 6) respectively.and where ȧρ is the mask of the scheme (4).The first two derivatives of (6) are given in (7) by using similar method as 17 . (4) (5) D (i) = 1, for i = 0, 0, for i � = 0,

Subdivision collocation method (SCM)
In this section, building upon the basic concepts of the subdivision scheme and its derivative calculations, we have developed a Subdivision Collocation Method (SCM) for solving the proposed problem.Additionally, this section includes the stability analysis and error estimations.

Numerical approximation
This section includes the detailed derivation of the Subdivision Collocation Method (SCM) for the proposed problem: Consider a partition of [0, L] that is equally divided by the knots Hence the approximate solution u(x, t) of (1) to exact solution U(x, t) based on the collocation approach can be expressed as is the basis function of (4) and C j (t) are the time dependent quantities which are to be evaluated.Now it is possible to simplify the approximation of the function u k i at the coordinates (x i , t k ) across the finite interval [x i , x i+1 ] by with P 2,i (x j ) = D ( ) and m = i = 0, 1, 2, . . ., N .To obtain the approximation of solution of values of P 2,i (x j ) and its derivative at the knots are needed.Since the values vanishes at the other knots.The approximations of the problem (1) at the time level at t j+1 + h can be considered by where successive time levels are k and k + 1 where k = 0, 1, 2, . . ., we will rearrange equations from (10)- (12) and the results will be displayed below after discretizing the first order time derivatives by forward difference scheme.
where t is the step size.Here Crank Nicloson approach is used to simplify the problem for θ = 0.5 , hence Eq. ( 13) takes the form where i = m = 0, 1, 2, 3, . . ., N .This above calculation leads to the generation of a linear system of order (N + 1) × (N + 9) with (N + 9) unknowns at each level of time.www.nature.com/scientificreports/now for i = m = 0 in Eq. ( 15), we have for i = m = 1 , Eq. ( 15), implies similarly continuing the sequence for i = m = N − 1 , Eq. ( 15), becomes for i = m = N , Eq. ( 15), implies from the above equations we got the system of equation and matrix form

Using the notation P
where After some simplification of W k+1 and W k with the help of (8), Eq. ( 20) can be transformed as Eq. ( 21), implies

Formation of consistent system
Since the order of the coefficient matrix of system ( 22) is (N + 1) × (N + 9) .After adding the given boundary conditions to the system (22), the order of the coefficient matrix increases to (N + 3) × (N + 9) .However, the system remains inconsistent.To attain a unique solution, six additional conditions are required.These conditions are constructed as follows: Since the degree of reproduction of subdivision scheme (4) is five, so we have established six-order boundary conditions to ensure a unique solution.For simplicity, we will only discuss the left end points, namely c −3 , c −2 , c −1 .The right end points c N+1 , c N+2 , c N+3 are handled in a similar manner.www.nature.com/scientificreports/ The values c −3 , c −2 , c −1 are determined by the quintic polynomial interpolating the points(u r , c r ) , 0 ≤ r ≤ 5 .i.e.The conditions at the left end are derived from where Since by (8), C(u r ) = c r for r = 1, 2, 3 and substituting u r by −u r in (24), we have By the equation given below we can get three left end conditions Using same pattern as left end conditions, we have following three right end conditions After obtaining these left and right end conditions the system of nonlinear equations of order (N + 9) × (N + 9) will be formed where where E is represented in (22), The first three rows of the matrix L c 0 are derived from all the requirements stated at the domain's left end from (25).And the fourth row of L c 0 is formed from (3) at C (0) = u 0 = φ 1 (t) .Similarly the last three rows of the matrix R c N are formed from all of the requirements that have been obtained at the right end of the domain from the (26).First row comes from (3) at C (N) = u L = φ 2 (t) .Hence The column vector M is defined as where C k+1 is given in (22).

Initial guess
Numerous iterative algorithms heavily depend on the initial guess, and it has a profound impact on the stability, speed, and reliability of the algorithm.When employing numerical and computational methods, the process of selecting a suitable initial guess is often intricate and problem-specific, and its significance cannot be underestimated.Selecting the appropriate initial guess, suited to the particular requirements of each scenario, is unquestionably a crucial step in tackling a wide range of mathematical and computational issues.
To initiate the iterative process (27), we develop an efficient method for computing the initial guess, which supports the convergence of iterative algorithms.The initial guess, denoted as is constructed from the initial and boundary conditions of the given problem as follows: At the initial time level, the approximate solution, as derived from ( 8), takes the form, where C ′ s s are unknown.For the determination of C s , the following conditions must be satisfied to ensure the correctness the initial approximation (30) as where x s = sh .As a result (N + 9) × (N + 9) matrix system where ψ 1 (x) , φ 1 (t) and φ 2 (t) and J c are given in (2), ( 3) and (28) respectively.Hence, the initial guess for the given problem will be obtained by solving Eq. (31).

Iterative algorithm
The iterative algorithm is defined in detail as follows: Step 1 Initial Solution: Calculate the initial solution, C 0 , by solving the system (31) using any numerical method.
The iterative algorithm is defined in the form of a flowchart given Fig. 1.The above equations can be written as where Using (37) and ( 38) in (36), we get From (39), it is concluded that the system is convergent for 4h 2 α t > 0 or α > 0 and Z > 0 the method.

Error estimation
The errors between the exact solutions U i and the numerical solutions u i are estimated at x i by calculating the absolute error , average error , L 2 , root mean square error, maximum relative error are defined as follows: The absolute error is denoted by Error ABS and estimated as The average error is denoted by Error AVG and is estimated as

Computational order
Computational order (C-order) of the subdivision collocation method presented with the following formula where E i−1 , E i and E i+1 are the errors at consecutive iterations, for grid size h.

Numerical examples, results and discussion
The numerical approach suggested in the preceding part is demonstrated here.In this section using the subdivision collocation algorithm to find solution of the different problems.Computations are carried out in Matlab.
The numerical solutions for each example are calculated as follows: we first computed an initial guess by solving (31), and then, using the obtained initial guess, the numerical solution are obtained through an iterative algorithm (32)-(33) for α > 0 , h > 0 and t > 0 .In Examples 1-5, if the advection time scale is less than the diffusion time scale the transport is advection controlled this corresponds P e >> 1 , and if P e << 1 then diffusion dominates.

Numerical examples
Example 1 We consider transport equation ( 1) with flow velocity β = 2 and diffusion coefficient |α| = |β|/P e , where P e = Pèclet number, in a channel length L = 2 of constant cross section and a bottom slop with physical conditions 1 initial condition: Example 2 We consider transport equation (1) taken from 1,23,24 with different values of flow velocity β , diffusion coefficient in channel length L of constant cross section and physical conditions, initial condition: and boundary conditions: φ 1 (t) = ae (bt) and φ 2 (t) = ae (bt−cL) .The exact solution is given by ae (bt−cx) , where .

Example 5
We consider transport equation (1) taken from 7,8 having L = 9 with constant cross section and bottom slope with the physical condition, ψ 1 (x) = exp(−(x − 1) 2 /α), and boundary conditions given below, and the exact solution is • The computational results of Example 1 are presented in Tables 1, 2, 3 for flow velocity β = 0.25 , t = 2h , L = 2 and N = 200 .Table 1 displays the comparison between the exact and approximate solutions for h = 0.01 , t = 2 , and P e = 2.5 , while also illustrating the accuracy and efficiency of the present algorithm through the absolute errors and CPU time = 0.337 s.Tables 2 and 3 display average error norms, computational order and CPU time for various choices of h with P e values of 2.5, 25, 250, and 2500 using the present method.They are compared with average errors reported in 1 .Both tables reveal that the present method yields good results for each value of P e .Figure 2 presents the exact and approximate solutions plot for P e = 2.5 , T = 1 , and h = 1 128 .From Fig. 2, we can observe that the numerical results of the present method closely align with the exact solutions.
• The computational results for Example 2 are presented in Tables 4, 5, 6 and 7 for various values of flow velocity β diffusion coefficient α , length L and P e .Table 4 displays the comparison between the exact and approximate solutions for β = −1 , L = 1 , h = 0.01 , t = 1 , and P e = 10 , while also illustrating the accuracy and efficiency of the present algorithm through the absolute errors and CPU time = 0.393 s for N .The Table 5 display the comparison between the exact and approximate solutions as well as comparison with existing methods 23,24 .The average errors for Example 2 obtained using the present method with parameters a = 1 , b = 0.1 , flow velocity β = −1 , time step t = 2h , L = 1 , and T = 1 are compared to existing methods 1 for various values of P e are presented in Tables 6 and 7.The data presented in the tables clearly indicates that Table 1.Comparison between the exact and approximate solutions of Example 1 using the Present SCM for P e = 2.5, h = 0.01, t = 2, t = 2h, L = 2.  www.nature.com/scientificreports/as values of P e increases, advection gains more control, resulting in smoother transport compared to the methods in 1,23,24 .A graph depicting the exact and approximate solutions for Example 2 at different values of P e is shown in Fig. 3.This comparison unequivocally establishes that the present method consistently provides significantly more accurate results.• The computational results for Example 3, with Dirichlet boundary conditions, are presented in Table 8.The results for flow velocity β = 1 , a diffusion coefficient of α = 1/1000 , a length of L = 2 , a spatial step size of h = 0.01 , a time interval of t = 1 , a time step size of t = 0.05h , and Pe = 1000 .The table compares the exact and approximate solutions, illustrating the accuracy and efficiency of the present algorithm through absolute errors.The CPU time for N = 100 is recorded at 0.393 s.The computation of absolute error results for Example 3 with Dirichlet boundary conditions, obtained using the present SCM, with the following parameters: flow velocity β = 1 , time step t = 0.05h , length L = 2 , and T = 1 , are compared to existing methods 2,6,23 for various values of P e are presented in Tables 9, 10. Figure 4 represents the exact and approximate solution plots for various values of P e = 0.5, 1.50, 2.0 , with h = 1/64 , and T = 1.

Grid points
The computation of absolute error L ∞ and L 2 for Example 3 with Neumann boundary conditions are pre- sented in Table 11.These results are obtained for a flow velocity of β = 0.1 , time step t = 0.1h , and a length of L = 2 .They are compared to the methods in 2,6 for different values of T. It is evident that the computational data and graphical representation of the current method closely align with the exact solution for different values of P e .
In both Dirichlet and Neumann boundary value problems, the transport is advection-controlled, rapid, and smooth, except in Table 9 where, for P e = 0.5 < 1 , the diffusion term dominates.
• The computational results of Example 4 are presented in Tables 12, 13, 14 for β = 1 and L = 1 and 2 .Table 12 displays the comparison between the exact and approximate solutions for h = 0.01 , t = 2 , L = 1 , N = 100 and P e = 10 , while also illustrating the accuracy and efficiency of the present algorithm through the absolute errors and CPU time = 0.339 s.The computation of absolute error results for Example 4 using the present method and a comparison of absolute errors with existing methods [1][2][3][4][5] are presented in Tables 13 and 14.These results are obtained for a diffusion coefficient of α = 0.1 , a flow velocity of β = 1 , a time step of t = 2h , a length of L = 1 , and a time duration of T = 2 , for different values of h = 0.1, 0.05, 0.01 .A graphical repre- sentation of the exact and approximate solutions for h = 0.01, L = 1, T = 2 is depicted in Fig. 5.
The numerical data and graphical representation closely match the exact solution.In Example 4, transportation is advection-controlled and more efficient than existing methods.
• The numerical solution for Example 5, obtained within a CPU time of 0.389 seconds using the present method with values of flow velocity β = 0.8 , diffusion coefficient α = 0.005 , t = 0.005 , h = 0.025 , T = 5 , and N = 100 , is presented in Table 15, and it is compared to existing methods 7,8 .The graph in Fig. 6 illustrates the exact and approximate solutions for Example 5 with h = 0.025 and T = 5 .Given these parameters, it is evident that the advection time scale is shorter than the diffusion time scale, indicating that transportation is advection-controlled and very fast.

Discussion on findings
In this research work, we took five examples from the literature and solve with current method and compare them with previously repointed methods other existing methods.In all examples, both computational and graphical results indicate that substances move rapidly rather than dispersing within the fluid.Additionally, it is observed that our current technique aligns well with the analytical solution compared to other methods.
• In Example 1, we compared our results with the numerical findings obtained by 1 .It is evident from the results presented in Tables 1, 2, 3 and 16 that the substances moved more swiftly compared to the results obtained by 1 , considering the constant flow velocity and diffusion coefficient.The substance concentrations are advected along the flow direction for the given parametric values, performing notably well within the specified interval 0 < β ≤ 0.25 for flow velocity and times 0 < t ≤ 2 .Employing our present method with smaller values of h allows for an increase in the constant flow velocity's magnitude, enabling swifter movement of substances within an extended channel.• Numerical results obtained from Example 2 are compared with those obtained from 1,23,24 , presented in Tables 4, 5, 6, 7 and 16.The computational data indicates a consistent fluid flow due to the swift movement of substances within the specified channel length, as evidenced in the results obtained from other existing methods.The current technique excels in addressing the problem within a time range of 0 < t ≤ 1 and a flow velocity range of −1 < β ≤ 1 , utilizing the same parameters as those employed in Example 2. By employing smaller h values and extending the channel length with a higher constant flow velocity, this research method achieves uniform flow and expedites the removal of substances present in the fluid within a minimal timeframe.• Tables 8, 9, 10 and 11 display the computational data of Example 3 along with the results obtained by the methods 2,6,23 .Table 16 represents the different types of errors obtained by present method.All these results demonstrate the accuracy of the present method, confirming a uniform flow of the fluid within the range of 0 < t ≤ 1 and 0 < β ≤ 1 for smaller values of h.Our research method yields favorable results within a short time period with a minimal velocity range for the same channel length.Moreover, it significantly improves the efficiency of moving substances within the fluid more rapidly and efficiently over longer distances with a high constant velocity.The method's efficiency increases notably with smaller values of h.• The numerical results of Example 4 are compared with those from [1][2][3][4][5] , presented in Tables 12, 13, 14.For the sake of comparison, the numerical results are obtained using the same grid points as those used in the existing methods.These results validate the swift transportation of substances and ensure a uniform flow of fluid within the specified channel length.Table 16 further validates the method's accuracy through various error evaluations for N = 100 .The method's accuracy remains consistent within the minimum ranges of 0 < t ≤ 2 and 0 < β ≤ 1 , remaining unaffected by variations.Moreover, this research method demonstrates improved performance for longer distances when employing increased constant flow velocity and smaller values of h.• The numerical results of Example 5, presented in Table 14, illustrate the uniform flow of the fluid based on the employed parameters.Additionally, Table 15 confirms the rapid movement of substances within the fluid compared to existing methods 7,8 .The accuracy of the current method for N = 100 is presented in Table 16, showcasing various error computations.The computational data demonstrates the efficient performance of the present method, confirming the rapid transport of substances within the limited ranges of 0 < t ≤ 5 and 0 < β ≤ 0.8 .Furthermore, increasing the channel length with smaller values of h and higher velocity does not significantly affect the accuracy of the present research method.

Conclusions
This article explores the application of the Subdivision Scheme-based Collocation Method (SCM) to numerically solve the one-dimensional advection-diffusion equation.A conventional finite difference approach is employed to discretize the time derivatives.The numerical results, presented in Tables 1-17 and Figs. 2, 3, 4, 5 and 6, demonstrate the reliability of the obtained outcomes.Specifically, the exact and approximate solutions for Examples 1 and 2 with absolute errors are represented in Tables 1 and 4. Average errors of Examples 1 and 2 presented in Tables 2, 3, 6, and 7 indicate that for larger values of P e , the current method is more efficient than the cubic C 1 spline-based method 1 .Table 5 represents the comparison between the exact and approximate solutions of Example 2 with the computational technique 23 , ICNTPS, EXPTPS, EULTPS methods 24 .Table 5 clearly demonstrates the accuracy of the current   9 and 10 for Example 3 show that for smaller to larger values of P e , our proposed method is more efficient than cubic B-Spline methods 1,2 and the computational technique 23 .www.nature.com/scientificreports/displays the computation of errors L ∞ and L 2 for different values of T, showing that SCM has significantly smaller errors than cubic B-Spline methods 2 and the methods of 6 .Table 12 for Example 4 represents the accuracy of the current method.Computations of Example 4 in Tables 13 and 14 for different values of h indicate the current method's accuracy for smaller values of h compared to the cubic B-spline [1][2][3][4] and exponential B-spline method 5 .
The data presented in Table 15 is the analysis of absolute errors of Example 5 and shows that the results obtained by SCM are very close to the exact solution when compared to the graph-theoretic polynomial collocation method 7 and the six-order compact finite difference method 8 .
Table 10.Comparison between the absolute errors of existing methods and the Present SCM of Example 3 at h = 1/64 , T = 1 and P e = 0.5, 1.5 and 2.0.Analysis of the data confirms the reliability, efficacy, and applicability of the current method.This technique can widely be used in fields such as fluid dynamics, heat transfer, and chemical engineering to understand and predict the behavior of different substances in various environments.In the future, we will explore this transport equation in two and three dimensions using the subdivision collocation method.

and estimation of errors Stability analysis Theorem 1
1, 2, 3, . .., The iterative algorithm (32) converges for 4h 2 α t > 0 or α > 0 and Z > 0.Proof The Von-Neumann stability approach provide us growth of error in as single Fourier mode.For this we use, Output: The Result at kth level when (33) satisfied Stop yes no yes no Figure 1.Flowchart of iterative algorithm.Vol.:(0123456789) Scientific Reports | (2024) 14:1712 | https://doi.org/10.1038/s41598-024-52059-7www.nature.com/scientificreports/Stability analysis The root mean square error is denoted by RMSE and is estimated as where N represents the total number of node points.The maximum relative error is denoted by MRE and is estimated as .

Table 2 .
Comparison between the average errors of existing methods and the Present SCM of Example 1 at P e = 2.5 and 25.
h for P

Table 3 .
Comparison between the average errors of existing methods and the Present SCM of Example 1 at P e = 250 and 2500.

Table 6 .
Comparison between the average errors of existing methods and the Present SCM of Example 2 at P e = 10 and 20.

Table 7 .
Comparison between the average errors of existing methods and the Present SCM of Example 2 at P e = 100 and 1000.
h for P

e = 100 for P e = 1000 SCM 1 C-order by SCM CPU time by SCM SCM 1 C-order by SCM CPU time by SCM
technique compared to the others.Exact and numerical solutions of Example 3 are listed in Table8which shows the accuracy of the current method.Tables

Table 8 .
Comparison between the exact and approximate solutions of Example 3 using the Present SCM for t = 0.05h, h = 0.01, t = 1, L = 2.

Table 9 .
Comparison between the absolute errors of existing methods and the Present SCM of Example 3 at h = 1/64 , T = 1 and P e = 1000, 10,000 and 20,000.

Table 11 .
Exact and approximate solutions of Example 3 at different P e numbers.Comparison between the L ∞ and L 2 errors of existing methods and the Present SCM of Example 3 at h = 1/10 , t = 0.1h.
T L

Table 16
shows the L ∞ , L 2 , Error AVG , RMSE, MRE of Examples 1-5 for N = 100.The comparison of the RMSE for all examples across different values of N is presented in Table 17, indicating that the RMSE remains relatively constant as N increases significantly.The graphical representations from Figs. 2, 3, 4, 5 and 6 depict that the transport model by SCM converges to the exact solution very quickly for different parameters in each case.

Table 12 .
Comparison between the exact and approximate solutions of Example 4 using the Present SCM for t = 2h, h = 0.01, t = 2, L = 1.